rm(list=ls())

setwd("MyDirectory")
source("Functions.R")
load("Data/RayonWeeks.RData")
ls()

summary(data)
names(data)
dim(data)
data <- streg.sort(data,unitvar="RID",timevar="WID")
head(data)

##########################
## Figure 0: Correlation matrix
##########################

names(data)
vars <- c("POP","FOREST","ELEVATION","SLOPE","GTD_SUICIDE","CHKDOM_NEAR","CHKINT_NEAR","DEPORTED","LANGUAGE","DIST_PIPES","REFUGEE_MIN","DIST_MIL","HOLIDAY","HOLIDAY_CHRI")
cor.mat <- cor(data[,vars],use="na.or.complete")
round(cor.mat,digits=1)
nog <- which(lower.tri(cor.mat)&(abs(cor.mat)>.4))
paste(rep(colnames(cor.mat),each=nrow(cor.mat)),colnames(cor.mat))[nog]

vars
varnames <- c("Population density","Forest","Elevation","Slope","Global suicide terrorism","Distance to internal checkpoint","Distance to international checkpoint","Deported in 1944","Percent Russian speaking","Distance to oil pipeline","Distance to nearest refugee camp","Distance to military base","Muslim holiday","ChRI holiday")
vars. <- vars[vars%in%names(data)]
varnames. <- varnames[vars%in%names(data)]
pdf("Output/Figure0_CorMat.pdf",width=5,height=4.5)
cormat(data=data,vars=vars.,labels=varnames.)
dev.off()


###########################
## Table 0: Summary stats
###########################

vars <- c("INS_ALL_I1","INS_ALL_N1","INS_ALL_I2","INS_ALL_N2","INS_ALL_I3","INS_ALL_N3","INS_ALL_I4","INS_ALL_N4","GOV_ALL","GOV_SEL","GOV_IND","INS_ALL_I1.b","INS_ALL_N1.b","INS_ALL_I2.b","INS_ALL_N2.b","INS_ALL_I3.b","INS_ALL_N3.b","INS_ALL_I4.b","INS_ALL_N4.b","GOV_ALL.b","GOV_SEL.b","GOV_IND.b","POP","FOREST","ELEVATION","SLOPE","GTD_SUICIDE","CHKDOM_NEAR","CHKINT_NEAR","DEPORTED","LANGUAGE","DIST_PIPES","REFUGEE_MIN","DIST_MIL","HOLIDAY","HOLIDAY_CHRI")
varnames <- c("Islamist violence (expansive)","Nationalist violence (expansive)","Islamist violence (intermediate)","Nationalist violence (intermediate)","Islamist violence (limited)","Nationalist violence (limited)","Islamist violence (target-based)","Nationalist violence (target-based)","Government Violence (all)","Government Violence (selective)","Government Violence (indiscriminate)","Islamist violence (expansive, binary)","Nationalist violence (expansive, binary)","Islamist violence (intermediate, binary)","Nationalist violence (intermediate, binary)","Islamist violence (limited, binary)","Nationalist violence (limited, binary)","Islamist violence (target-based, binary)","Nationalist violence (target-based, binary)","Government Violence (all, binary)","Government Violence (selective, binary)","Government Violence (indiscriminate, binary)","Population density","Forest","Elevation","Slope","Global suicide terrorism","Distance to internal checkpoint","Distance to international checkpoint","Deported in 1944","Percent Russian speaking","Distance to oil pipeline","Distance to nearest refugee camp","Distance to military base","Muslim holiday","ChRI holiday")

summarystat <- function(x){
  range <- paste("[",round(min(x,na.rm=T),3),", ",round(max(x,na.rm=T),3),"]",sep="")
  medz <- round(median(x,na.rm=T),3)
  meanz <- round(mean(x,na.rm=T),3)
  sdz <- round(sd(x,na.rm=T),3)
  return(c(range,medz,meanz,sdz))
}

summat <- t(apply(data[,vars],2,summarystat))
row.names(summat) <- varnames
colnames(summat) <- c("Range","Median","Mean","Std. Dev.")
xtable(summat)
write.table(summat,file="Output/Table0_SumStats.txt")
write.table(print(xtable(summat)),file="Output/Table0_SumStats.tex")



####################
## Figure 2: Time plots
####################

length(c(1:31,1:28,1:31,1:30,1:31,1:30))/7

# Create time series
ts.data <-aggregate(data[,c("GOV_ALL","GOV_ALL_F","GOV_ALL_L","GOV_ALL_O","INS_ALL1","INS_ALL_I1","INS_ALL_N1","INS_ALL_O1","INS_ALL2","INS_ALL_I2","INS_ALL_N2","INS_ALL_O2","INS_ALL3","INS_ALL_I3","INS_ALL_N3","INS_ALL_O3","INS_ALL4","INS_ALL_I4","INS_ALL_N4","INS_ALL_O4")],by=list(WID=data$WID),FUN=sum)

## Plot
defz <- c("Expanded definition","Intermediate definition","Limited definition","Target-based")
i <- 1
for(i in 1:4){
# Absolute numbers plot
xx1 <- ts.data[,paste("INS_ALL_N",i,sep="")]
ts1 <- ts(xx1, start=c(2000,26), frequency=52)
ts1 <- stl(ts1, s.window="periodic")
xx2 <- ts.data[,paste("INS_ALL_I",i,sep="")]
ts2 <- ts(xx2, start=c(2000,26), frequency=52)
ts2 <- stl(ts2, s.window="periodic")
jpeg(width=9,height=1.5,units="in",file=paste("Output/Figure2_Timeline_",i,"a.jpg",sep=""),res=300)
par(mar=c(2,4,.5,.5))
plot(ts.data$WID,xx1,main=NA,xlab="Time",ylab="Events",xaxt="n",cex.axis=.7,cex=.3,pch=3,col=NA,las=1,bty="l")
points(ts.data$WID,xx1,cex=.3,pch=4,col="red")
points(ts.data$WID,xx2,cex=.3,pch=4,col="darkgreen")
lines(c(ts1[[1]][,2]),col="black",lwd=2)
lines(c(ts1[[1]][,2]),col="red")
lines(c(ts2[[1]][,2]),col="black",lwd=2)
lines(c(ts2[[1]][,2]),col="green")
axis(1,at=((0:11)*52+26),labels=2001:2012,cex.axis=.7)
legend(x="topright",pch=c(3,4),col=c("red","darkgreen"),legend=c("Nationalist","Islamist"),bg="white",cex=.65,title=defz[i],ncol=2,bty="n")
dev.off()

## Prortional, shaded area
xx <- ts.data[,paste("INS_ALL_I",i,sep="")]/(ts.data[,paste("INS_ALL_I",i,sep="")]+ts.data[,paste("INS_ALL_N",i,sep="")])
xx[is.na(xx)] <- 0
ts1 <- ts(xx[1:621], start=c(2000,26), frequency=52)
ts1 <- stl(ts1, s.window="periodic")
jpeg(width=9,height=.8,units="in",file=paste("Output/Figure2_Timeline_",i,"b.jpg",sep=""),res=300)
par(mar=c(0.4,4,0.2,.5))
plot(ts.data$WID,xx,main=NA,xlab="Time",ylab="Proportion",yaxt="n",xaxt="n",cex.axis=.7,cex=.3,pch=4,col=NA,las=1,bty="l",ylim=c(0,.5))
#axis(1,at=(((-1):11)*52+26),labels=c(NA,2001:2012),cex.axis=.7)
axis(2,cex.axis=.7,las=1)
#polygon(c(ts.data$WID[c(1,1:621,621)]),c(1,ts1[[1]][,2],1),col="pink",border=NA)
#segments(x0=-1,x1=nrow(ts.data),y0=(0:5)*.1,y1=(0:5)*.1,col="grey",lty="dotted")
polygon(c(ts.data$WID[c(1:621,621:1)]),c(ts1[[1]][,2],rep(0,621)),col="darkgreen",border=NA)
par(xpd=NA)
segments(x0=(((0):11)*52+26),x1=(((0):11)*52+26),y0=1.2,y1=.45)
dev.off()
}



######################
## Table 1: by type
######################

sum.mat <- data.frame(Islamist=rep(NA,4),Nationalist=rep(NA,4),Other=rep(NA,4),Total=rep(NA,4))
row.names(sum.mat) <- c("Expansive","Intermediate","Limited","Target-based")


i <- 1
for(i in 1:4){
  sum.mat[i,1] <- paste(sum(ts.data[,paste("INS_ALL_I",i,sep="")])," (",round((100*sum(ts.data[,paste("INS_ALL_I",i,sep="")])/(sum(ts.data[,paste("INS_ALL_I",i,sep="")])+sum(ts.data[,paste("INS_ALL_N",i,sep="")])+sum(ts.data[,paste("INS_ALL_O",i,sep="")]))),1),"%)",sep="")
  sum.mat[i,2] <- paste(sum(ts.data[,paste("INS_ALL_N",i,sep="")])," (",round((100*sum(ts.data[,paste("INS_ALL_N",i,sep="")])/(sum(ts.data[,paste("INS_ALL_I",i,sep="")])+sum(ts.data[,paste("INS_ALL_N",i,sep="")])+sum(ts.data[,paste("INS_ALL_O",i,sep="")]))),1),"%)",sep="")
  sum.mat[i,3] <- paste(sum(ts.data[,paste("INS_ALL_O",i,sep="")])," (",round((100*sum(ts.data[,paste("INS_ALL_O",i,sep="")])/(sum(ts.data[,paste("INS_ALL_I",i,sep="")])+sum(ts.data[,paste("INS_ALL_N",i,sep="")])+sum(ts.data[,paste("INS_ALL_O",i,sep="")]))),1),"%)",sep="")
  sum.mat[i,4] <- paste(sum(ts.data[,paste("INS_ALL_I",i,sep="")])+sum(ts.data[,paste("INS_ALL_N",i,sep="")])+sum(ts.data[,paste("INS_ALL_O",i,sep="")])," (100%)",sep="")
}
sum.mat
write.table(sum.mat,file="Output/Table1_ByType.txt")
write.table(print(xtable(sum.mat)),file="Output/Table1_ByType.tex")


######################
## Figure 3a: by oblast
######################

i <- 2
obl.mat <- aggregate(data[,c(paste("INS_ALL_N",i,sep=""),paste("INS_ALL_I",i,sep=""))],by=list(GROUP=data$OBLAST_NM2),FUN=sum)
obl.mat[,2] <- obl.mat[,2]/sum(data[,paste("INS_ALL_N",i,sep="")])*100
obl.mat[,3] <- obl.mat[,3]/sum(data[,paste("INS_ALL_I",i,sep="")])*100
obl.mat
oblasts <- obl.mat[obl.mat[,2]>0|obl.mat[,3]>0,]
oblasts[oblasts[,1]=="OTHER",-1] <- apply(oblasts[!oblasts[,1]=="OTHER",-1],2,function(x){100-sum(x)})
oblasts <- oblasts[order(oblasts[,2],decreasing=T),]
write.table(oblasts,file="Output/Table2_ByOblast.txt")
write.table(print(xtable(oblasts)),file="Output/Table2_ByOblast.tex")

## Plot
cols <- c("white","grey95","grey95","grey15","grey15","grey10","grey10")
pdf(file="Output/Figure3a_oblast.pdf",height=3,width=4)
par(mar=c(2,4,.5,.5))
plot(c(0,1),c(0,100),col=NA,xaxt="n",xlab="",ylab="Percent of total",bty="l")
axis(1,at=c(.25,.75),label=c("Nationalist","Islamist"))
prop.plot2(obl.mat[,2],toupper(obl.mat[,1]),n=7,x0 = .05, x1 =.45,lab.print=T,cols=cols,other=F)
prop.plot2(obl.mat[,3],toupper(obl.mat[,1]),n=7,x0 = .55, x1 =.95,lab.print=T,cols=cols,other=F)
dev.off()


######################
## Figure 3b: by ethnicity
######################

i <- 2
eth.mat <- aggregate(data[,c(paste("INS_ALL_N",i,sep=""),paste("INS_ALL_I",i,sep=""))],by=list(GROUP=data$GROUP1),FUN=sum)
eth.mat[,2] <- eth.mat[,2]/sum(data[,paste("INS_ALL_N",i,sep="")])*100
eth.mat[,3] <- eth.mat[,3]/sum(data[,paste("INS_ALL_I",i,sep="")])*100
eth.mat
ethnicities <- eth.mat[eth.mat[,2]>1|eth.mat[,3]>1,]
ethnicities[ethnicities[,1]=="OTHER",-1] <- apply(ethnicities[!ethnicities[,1]=="OTHER",-1],2,function(x){100-sum(x)})
ethnicities <- ethnicities[order(ethnicities[,2],decreasing=T),]
write.table(ethnicities,file="Output/Table3_ByEth.txt")
write.table(print(xtable(ethnicities)),file="Output/Table3_ByEth.tex")


## Plot
cols <- c("white","grey95","grey95","grey95","grey15","grey15","grey15","grey15","grey25","grey10","black")
pdf(file="Output/Figure3b_ethnicity.pdf",height=3,width=4)
par(mar=c(2,4,.5,.5))
plot(c(0,1),c(0,100),col=NA,xaxt="n",xlab="",ylab="Percent of total",bty="l")
axis(1,at=c(.25,.75),label=c("Nationalist","Islamist"))
prop.plot2(eth.mat[,2],eth.mat[,1],n=8,x0 = .05, x1 =.45,lab.print=T,cols=cols,other=T)
prop.plot2(eth.mat[,3],eth.mat[,1],n=10,x0 = .55, x1 =.95,lab.print=T,cols=cols,other=T)
dev.off()




######################
## Table 0: proportion by year
######################
names(data)

i <- 1
years <- 2000:2012
prop.t <- as.data.frame(matrix(NA,length(years),5))
names(prop.t) <- c("Year","Expansive","Intermediate","Limited","Target-based")
prop.t[,1] <- years

for(i in 1:length(years)){
  prop.t[i,2] <- sum(data[data$YEAR==years[i],"INS_ALL_I1"],na.rm=T)/sum(data[data$YEAR==years[i],"INS_ALL"],na.rm=T)*100
  prop.t[i,3] <- sum(data[data$YEAR==years[i],"INS_ALL_I2"],na.rm=T)/sum(data[data$YEAR==years[i],"INS_ALL"],na.rm=T)*100
  prop.t[i,4] <- sum(data[data$YEAR==years[i],"INS_ALL_I3"],na.rm=T)/sum(data[data$YEAR==years[i],"INS_ALL"],na.rm=T)*100
  prop.t[i,5] <- sum(data[data$YEAR==years[i],"INS_ALL_I4"],na.rm=T)/sum(data[data$YEAR==years[i],"INS_ALL"],na.rm=T)*100
}
xtable(prop.t,digits=2)

write.table(prop.t,file="Output/Table0_ByYear.txt")
write.table(print(xtable(prop.t)),file="Output/Table0_ByYear.tex")



######################
## Table 2: Distance from attack
######################


# Load village-level data
rm(data)
load("Data/CaucasusData.RData")
names(data)
dim(data)

dist.mat <- as.data.frame(matrix(NA,4,3))
names(dist.mat) <- c("Islamist","Nationalist","KS Test")
rownames(dist.mat) <- c("Expanded","Intermediate","Limited","Target-based")

## CI
dist.mat <- as.data.frame(matrix(NA,4,4))
names(dist.mat) <- c("Islamist","Islamist CI","Nationalist","Nationalist CI")
rownames(dist.mat) <- c("Expanded","Intermediate","Limited","Target-based")

# 0=Limited, 1=Intermediate, 2=Expansive, 3=Target-based
i <- 2
sum(data[,paste("INS_ALL_I",i,sep="")],na.rm=T)

for(i in 1:4){
  dist.mat[i,1] <- mean(data[data[,paste("INS_ALL_I",c(2,1,0,3)[i],sep="")]>0,paste("L_W_INS_ALL_I",c(2,1,0,3)[i],sep="")],na.rm=T) 
  x <- data[data[,paste("INS_ALL_I",c(2,1,0,3)[i],sep="")]>0,paste("L_W_INS_ALL_I",c(2,1,0,3)[i],sep="")]
  bootz <- c()
  for(k in 1:1000){
  	bootz[k] <- mean(x[sample(1:length(x),length(x),replace=T)],na.rm=T)}
dist.mat[i,2] <- paste("(",round(quantile(bootz,.025),2),", ",round(quantile(bootz,.975),2),")",sep="")

  dist.mat[i,3] <- mean(data[data[,paste("INS_ALL_N",c(2,1,0,3)[i],sep="")]>0,paste("L_W_INS_ALL_N",c(2,1,0,3)[i],sep="")],na.rm=T)
    x <- data[data[,paste("INS_ALL_N",c(2,1,0,3)[i],sep="")]>0,paste("L_W_INS_ALL_N",c(2,1,0,3)[i],sep="")]
  bootz <- c()
  for(k in 1:1000){
  	bootz[k] <- mean(x[sample(1:length(x),length(x),replace=T)],na.rm=T)}
dist.mat[i,4] <- paste("(",round(quantile(bootz,.025),2),", ",round(quantile(bootz,.975),2),")",sep="")
  
#  t. <- ks.test(data[data[,paste("INS_ALL_I",c(2,1,0,3)[i],sep="")]>0,paste("L_W_INS_ALL_I",c(2,1,0,3)[i],sep="")],data[data[,paste("INS_ALL_N",c(2,1,0,3)[i],sep="")]>0,paste("L_W_INS_ALL_N",c(2,1,0,3)[i],sep="")])
#  dist.mat[i,3] <-  paste(round(t.$statistic,2),ifelse(t.$p.value<.001,"***",ifelse(t.$p.value<.01,"**",ifelse(t.$p.value<.05,"*",ifelse(t.$p.value<.1,"'","")))),sep="")
}
dist.mat

write.table(dist.mat,file="Output/Table2_DistanceToNearest.txt")
write.table(print(xtable(dist.mat)),file="Output/Table2_DistanceToNearest.tex")





######################
## Table 3: GTD
######################

rm(data)
load("Data/RayonWeeks.RData")

gtd.mat <- as.data.frame(matrix(NA,4,5))
names(gtd.mat) <- c("Islamist","t-test (I)","Nationalist","t-test (N)","t-test (I vs. N)")
rownames(gtd.mat) <- c("Expanded","Intermediate","Limited","Target-based")
i <-1

## CI
gtd.mat <- as.data.frame(matrix(NA,4,4))
names(gtd.mat) <- c("Islamist","Islamist CI","Nationalist","Nationalist CI")
rownames(gtd.mat) <- c("Expanded","Intermediate","Limited","Target-based")


for(i in 1:4){
  gtd.mat[i,1] <- 100*(mean(data[data$GTD_SUICIDE>mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_I",i,sep="")])/mean(data[data$GTD_SUICIDE<=mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_I",i,sep="")])-1)
  #  t. <- t.test(data[data$GTD_SUICIDE<=mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_I",i,sep="")],data[data$GTD_SUICIDE>mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_I",i,sep="")])
#  gtd.mat[i,2] <- paste("(",round(t.$statistic,2),")",ifelse(t.$p.value<.001,"***",ifelse(t.$p.value<.01,"**",ifelse(t.$p.value<.05,"*",ifelse(t.$p.value<.1,"'","")))),sep="")
  	x1 <- data[data$GTD_SUICIDE>mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_I",i,sep="")]
  	x0 <- data[data$GTD_SUICIDE<=mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_I",i,sep="")]
  bootz <- c()
  for(k in 1:1000){
  	bootz[k] <- 100*(mean(x1[sample(1:length(x1),length(x1),replace=T)],na.rm=T)/mean(x0[sample(1:length(x0),length(x0),replace=T)],na.rm=T)-1)}
gtd.mat[i,2] <- paste("(",round(quantile(bootz,.025),2),", ",round(quantile(bootz,.975),2),")",sep="")
gtd.mat[i,3] <- 100*(mean(data[data$GTD_SUICIDE>mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_N",i,sep="")])/mean(data[data$GTD_SUICIDE<=mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_N",i,sep="")])-1)
#  t. <- t.test(data[data$GTD_SUICIDE<=mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_N",i,sep="")],data[data$GTD_SUICIDE>mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_N",i,sep="")])
#  gtd.mat[i,4] <- paste("(",round(t.$statistic,2),")",ifelse(t.$p.value<.001,"***",ifelse(t.$p.value<.01,"**",ifelse(t.$p.value<.05,"*",ifelse(t.$p.value<.1,"'","")))),sep="")
  	x1 <- data[data$GTD_SUICIDE>mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_N",i,sep="")]
  	x0 <- data[data$GTD_SUICIDE<=mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_N",i,sep="")]
  bootz <- c()
  for(k in 1:1000){
  	bootz[k] <- 100*(mean(x1[sample(1:length(x1),length(x1),replace=T)],na.rm=T)/mean(x0[sample(1:length(x0),length(x0),replace=T)],na.rm=T)-1)}
gtd.mat[i,4] <- paste("(",round(quantile(bootz,.025),2),", ",round(quantile(bootz,.975),2),")",sep="")
#  t. <- t.test(data[data$GTD_SUICIDE>mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_I",i,sep="")],data[data$GTD_SUICIDE>mean(data$GTD_SUICIDE,na.rm=T),paste("INS_ALL_N",i,sep="")])
#  gtd.mat[i,5] <- paste("(",round(t.$statistic,2),")",ifelse(t.$p.value<.001,"***",ifelse(t.$p.value<.01,"**",ifelse(t.$p.value<.05,"*",ifelse(t.$p.value<.1,"'","")))),sep="")
}

gtd.mat
write.table(gtd.mat,file="Output/Table3_GTD.txt")
write.table(print(xtable(gtd.mat)),file="Output/Table3_GTD.tex")



######################
## Figure 4: Muslim holidays
######################

names(data)
holi.mat <- as.data.frame(matrix(NA,4,5))
names(holi.mat) <- c("Islamist","t-test (I)","Nationalist","t-test (N)","t-test (I vs. N)")
rownames(holi.mat) <- c("Expanded","Intermediate","Limited","Target-based")
i <-1


## CI
holi.mat <- as.data.frame(matrix(NA,4,4))
names(holi.mat) <- c("Islamist","Islamist CI","Nationalist","Nationalist CI")
rownames(holi.mat) <- c("Expanded","Intermediate","Limited","Target-based")


for(i in 1:4){
  holi.mat[i,1] <- 100*(mean(data[data$HOLIDAY==1,paste("INS_ALL_I",i,sep="")])/mean(data[data$HOLIDAY==0,paste("INS_ALL_I",i,sep="")])-1)
#  t. <- t.test(data[data$HOLIDAY==1,paste("INS_ALL_I",i,sep="")],data[data$HOLIDAY==0,paste("INS_ALL_I",i,sep="")])
#  holi.mat[i,2] <- paste("(",round(t.$statistic,2),")",ifelse(t.$p.value<.001,"***",ifelse(t.$p.value<.01,"**",ifelse(t.$p.value<.05,"*",ifelse(t.$p.value<.1,"'","")))),sep="")
  	x1 <- data[data$HOLIDAY==1,paste("INS_ALL_I",i,sep="")]
  	x0 <- data[data$HOLIDAY==0,paste("INS_ALL_I",i,sep="")]
  bootz <- c()
  for(k in 1:1000){
  	bootz[k] <- 100*(mean(x1[sample(1:length(x1),length(x1),replace=T)],na.rm=T)/mean(x0[sample(1:length(x0),length(x0),replace=T)],na.rm=T)-1)}
holi.mat[i,2] <- paste("(",round(quantile(bootz,.025),2),", ",round(quantile(bootz,.975),2),")",sep="")
  holi.mat[i,3] <- 100*(mean(data[data$HOLIDAY==1,paste("INS_ALL_N",i,sep="")])/mean(data[data$HOLIDAY==0,paste("INS_ALL_N",i,sep="")])-1)
#  t. <- t.test(data[data$HOLIDAY==1,paste("INS_ALL_N",i,sep="")],data[data$HOLIDAY==0,paste("INS_ALL_N",i,sep="")])
#  holi.mat[i,4] <- paste("(",round(t.$statistic,2),")",ifelse(t.$p.value<.001,"***",ifelse(t.$p.value<.01,"**",ifelse(t.$p.value<.05,"*",ifelse(t.$p.value<.1,"'","")))),sep="")
  	x1 <- data[data$HOLIDAY==1,paste("INS_ALL_N",i,sep="")]
  	x0 <- data[data$HOLIDAY==0,paste("INS_ALL_N",i,sep="")]
  bootz <- c()
  for(k in 1:1000){
  	bootz[k] <- 100*(mean(x1[sample(1:length(x1),length(x1),replace=T)],na.rm=T)/mean(x0[sample(1:length(x0),length(x0),replace=T)],na.rm=T)-1)}
holi.mat[i,4] <- paste("(",round(quantile(bootz,.025),2),", ",round(quantile(bootz,.975),2),")",sep="")
#  t. <- t.test(data[data$HOLIDAY==1,paste("INS_ALL_I",i,sep="")],data[data$HOLIDAY==1,paste("INS_ALL_N",i,sep="")])
#  holi.mat[i,5] <- paste("(",round(t.$statistic,2),")",ifelse(t.$p.value<.001,"***",ifelse(t.$p.value<.01,"**",ifelse(t.$p.value<.05,"*",ifelse(t.$p.value<.1,"'","")))),sep="")
}

holi.mat
write.table(holi.mat,file="Output/Table4a_HOLIDAY.txt")
write.table(print(xtable(holi.mat)),file="Output/Table4a_HOLIDAY.tex")






######################
## Table 5: selective violence
######################

summary(data)

sel.mat <- as.data.frame(matrix(NA,4,4))
names(sel.mat) <- c("Islamist","Islamist CI","Nationalist","Nationalist CI")
rownames(sel.mat) <- c("Expanded","Intermediate","Limited","Target-based")
j <- 2
for(j in 1:4){
  
  predata <- prematch.data(treatment="SI",lag=12,def=j,type="I",bin=1,data=data)
  predata2 <- predata[,c("Treat","Post","Pre")]
  names(predata2)[2:3] <- c("Post_I","Pre_I")
  predata <- prematch.data(treatment="SI",lag=12,def=j,type="N",bin=1,data=data)
  predata2$Post_N <- predata$Post
  predata2$Pre_N <- predata$Pre
  
  # Bootstrap CI
  boots <- data.frame(DD_I=rep(NA,1000),DD_N=rep(NA,1000))
  for(i in 1:1000){
    sampT <- sample(which(predata2$Treat==1),sum(predata2$Treat==1),replace=T)
    sampC <- sample(which(predata2$Treat==0),sum(predata2$Treat==0),replace=T)
    boots[i,"DD_I"] <- 100*((mean(predata2$Post_I[sampT],na.rm=T)/mean(predata2$Post_I[sampC],na.rm=T))-1)
    boots[i,"DD_N"] <- 100*((mean(predata2$Post_N[sampT],na.rm=T)/mean(predata2$Post_N[sampC],na.rm=T))-1)
  }
  sel.mat[j,1] <- 100*((mean(predata2$Post_I[predata2$Treat==1],na.rm=T)/mean(predata2$Post_I[predata2$Treat==0],na.rm=T))-1)
  sel.mat[j,3] <- 100*((mean(predata2$Post_N[predata2$Treat==1],na.rm=T)/mean(predata2$Post_N[predata2$Treat==0],na.rm=T))-1)
  sel.mat[j,2] <- paste("(",round(quantile(boots[,"DD_I"],.025),2),", ",round(quantile(boots[,"DD_I"],.975),2),")",sep="")
  sel.mat[j,4] <- paste("(",round(quantile(boots[,"DD_N"],.025),2),", ",round(quantile(boots[,"DD_N"],.975),2),")",sep="")
}

sel.mat
write.table(sel.mat,file="Output/Table5_SEL.txt")
write.table(print(xtable(sel.mat)),file="Output/Table5_SEL.tex")




##########################
## Table 6
##########################


# Islamist violence
mod_i2a <- glm(INS_ALL_I2.b ~ L_INS_ALL_I2.b + L_Wb_INS_ALL_I2 + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI + POP + SLOPE + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_i2a)

mod_i2b <- glm(INS_ALL_I2.b ~ L_INS_ALL_I2.b + L_Wb_INS_ALL_I2 + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI + POP + ELEVATION + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_i2b)

mod_i2c <- glm(INS_ALL_I2.b ~ L_INS_ALL_I2.b + L_Wb_INS_ALL_I2  + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI + POP + LANGUAGE + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_i2c)

mod_i2d <- glm(INS_ALL_I2.b ~ L_INS_ALL_I2.b + L_Wb_INS_ALL_I2 + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI + POP + CHKINT_NEAR + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_i2d)

mod_i2e <- glm(INS_ALL_I2.b ~ L_INS_ALL_I2.b + L_Wb_INS_ALL_I2  + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI  + POP + DIST_PIPES + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_i2e)

mod_i2f <- glm(INS_ALL_I2.b ~ L_INS_ALL_I2.b + L_Wb_INS_ALL_I2  + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI  + POP + FOREST + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_i2f)

mod_i2g <- glm(INS_ALL_I2.b ~ L_INS_ALL_I2.b + L_Wb_INS_ALL_I2 + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI + POP + CHKINT_NEAR + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_i2g)


# Nationalist violence
mod_n2a <- glm(INS_ALL_N2.b ~ L_INS_ALL_N2.b + L_Wb_INS_ALL_N2  + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI  + POP + SLOPE + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_n2a)

mod_n2b <- glm(INS_ALL_N2.b ~ L_INS_ALL_N2.b + L_Wb_INS_ALL_N2  + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI  + POP + ELEVATION + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_n2b)

mod_n2c <- glm(INS_ALL_N2.b ~ L_INS_ALL_N2.b + L_Wb_INS_ALL_N2  + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI  + POP + LANGUAGE + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_n2c)

mod_n2d <- glm(INS_ALL_N2.b ~ L_INS_ALL_N2.b + L_Wb_INS_ALL_N2 + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI + POP + CHKINT_NEAR + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_n2d)

mod_n2e <- glm(INS_ALL_N2.b ~ L_INS_ALL_N2.b + L_Wb_INS_ALL_N2  + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI  + POP + DIST_PIPES + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_n2e)

mod_n2f <- glm(INS_ALL_N2.b ~ L_INS_ALL_N2.b + L_Wb_INS_ALL_N2 + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI + POP + FOREST + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_n2f)

mod_n2g <- glm(INS_ALL_N2.b ~ L_INS_ALL_N2.b + L_Wb_INS_ALL_N2 + GTD_SUICIDE + HOLIDAY + HOLIDAY_CHRI + POP + CHKINT_NEAR + DIST_MIL + REFUGEE_MIN + DEPORTED, data=data, family="binomial")
summary(mod_n2g)

## CF plots

vals <- quantile(data$GTD_SUICIDE,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2d,var="GTD_SUICIDE",vals=vals)
pred.c <- predman2(mod_n2d,var="GTD_SUICIDE",vals=vals)
jpeg(file="Output/Table6_GTD_SUICIDE.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-25,275),.8)
dev.off()

#vals <- quantile(data$HOLIDAY,c(.01,.99),na.rm=T)
vals <- c(0,1)
pred.b <- predman2(mod_i2f,var="HOLIDAY_CHRI",vals=vals)
pred.c <- predman2(mod_n2f,var="HOLIDAY_CHRI",vals=vals)
jpeg(file="Output/Table6_HOLIDAY_CHRI.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-50,50),.8)
dev.off()

#vals <- quantile(data$HOLIDAY,c(.01,.99),na.rm=T)
vals <- c(0,1)
pred.b <- predman2(mod_i2d,var="HOLIDAY",vals=vals)
pred.c <- predman2(mod_n2d,var="HOLIDAY",vals=vals)
jpeg(file="Output/Table6_HOLIDAY.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-50,50),.8)
dev.off()

vals1 <- quantile(data$L_Wb_INS_ALL_I2,c(.01,.99),na.rm=T)
vals2 <- quantile(data$L_Wb_INS_ALL_N2,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2d,var="L_Wb_INS_ALL_I2",vals=vals1)
pred.c <- predman2(mod_n2d,var="L_Wb_INS_ALL_N2",vals=vals2)
jpeg(file="Output/Table6_L_Wb_INS_ALL.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-25,175),.8)
dev.off()

vals <- quantile(data$SLOPE,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2a,var="SLOPE",vals=vals)
pred.c <- predman2(mod_n2a,var="SLOPE",vals=vals)
jpeg(file="Output/Table6_SLOPE.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-25,100),.8)
dev.off()

vals <- quantile(data$ELEVATION,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2b,var="ELEVATION",vals=vals)
pred.c <- predman2(mod_n2b,var="ELEVATION",vals=vals)
jpeg(file="Output/Table6_ELEVATION.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-25,100),.8)
dev.off()

vals <- quantile(data$LANGUAGE,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2c,var="LANGUAGE",vals=vals)
pred.c <- predman2(mod_n2c,var="LANGUAGE",vals=vals)
jpeg(file="Output/Table6_LANGUAGE.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-75,25),.8)
dev.off()

vals <- quantile(data$CHKINT_NEAR,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2d,var="CHKINT_NEAR",vals=vals)
pred.c <- predman2(mod_n2d,var="CHKINT_NEAR",vals=vals)
jpeg(file="Output/Table6_CHKINT_NEAR.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-100,25),.8)
dev.off()

vals <- quantile(data$DIST_PIPES,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2e,var="DIST_PIPES",vals=vals)
pred.c <- predman2(mod_n2e,var="DIST_PIPES",vals=vals)
jpeg(file="Output/Table6_DIST_PIPES.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-50,50),.8)
dev.off()

vals <- quantile(data$POP,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2d,var="POP",vals=vals)
pred.c <- predman2(mod_n2d,var="POP",vals=vals)
jpeg(file="Output/Table6_POP.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-25,200),.8)
dev.off()

vals <- quantile(data$DIST_MIL,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2d,var="DIST_MIL",vals=vals)
pred.c <- predman2(mod_n2d,var="DIST_MIL",vals=vals)
jpeg(file="Output/Table6_DIST_MIL.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-50,50),.8)
dev.off()

vals <- quantile(data$REFUGEE_MIN,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2d,var="REFUGEE_MIN",vals=vals)
pred.c <- predman2(mod_n2d,var="REFUGEE_MIN",vals=vals)
jpeg(file="Output/Table6_REFUGEE_MIN.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-100,25),.8)
dev.off()

vals <- quantile(data$DEPORTED,c(.01,.99),na.rm=T)
vals <- c(0,1)
pred.b <- predman2(mod_i2d,var="DEPORTED",vals=vals)
pred.c <- predman2(mod_n2d,var="DEPORTED",vals=vals)
jpeg(file="Output/Table6_DEPORTED.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-25,300),.8)
dev.off()

vals <- quantile(data$FOREST,c(.01,.99),na.rm=T)
pred.b <- predman2(mod_i2f,var="FOREST",vals=vals)
pred.c <- predman2(mod_n2f,var="FOREST",vals=vals)
jpeg(file="Output/Table6_FOREST.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot2(pred.b,pred.c,c(-25,100),.8)
dev.off()







########################################################
########################################################
########################################################
########################################################
## Counterinsurgency
########################################################
########################################################
########################################################
########################################################



##################################
## Diff-in-Diff setup
##################################

names(data)
data$FC <- as.numeric(data$GOV_ALL_F>0)
data$GOV_ALL_F.b <- as.numeric(data$GOV_ALL_F>0)
data$GOV_ALL_L.b <- as.numeric(data$GOV_ALL_L>0)
data$GOV_IND.b <- as.numeric(data$GOV_IND>0)
data$GOV_SEL.b <- as.numeric(data$GOV_SEL>0)


######################################
## Selective (T) vs. Indiscriminate (C)
######################################

rightside <- c("Pre","PreG","POP","REFUGEE_MIN","DIST_MIL","DIST_PIPES","DEPORTED","ELEVATION","SLOPE","LANGUAGE","FOREST","CHKINT_NEAR","GTD_SUICIDE","HOLIDAY","LONG","LAT","YEAR","MONTH","WID")
#rightside <- c("Pre","PreG","ELEVATION","FOREST","POP","LANGUAGE","REFUGEE_MIN","DIST_MIL","DEPORTED","YEAR","MONTH","LONG","LAT","CHKINT_NEAR","GTD_SUICIDE","HOLIDAY","HOLIDAY_CHRI")
form <- as.formula(paste("Treat ~", paste(rightside,collapse="+")))
exacts <- c("YEAR","MONTH","HOLIDAY")



###########################
## CEM match
###########################

rightside <- c("Pre","PreG","POP","REFUGEE_MIN","DIST_MIL","DIST_PIPES","DEPORTED","ELEVATION","SLOPE","LANGUAGE","FOREST","CHKINT_NEAR","GTD_SUICIDE","HOLIDAY","LONG","LAT","YEAR","MONTH","WID")

predata.n <- prematch.data(treatment="SI",lag=12,def=2,type="N",bin=1,data=data)
X <- predata.n[,c("Post","Treat",rightside)]
X <- na.omit(X)
mat <- cem(treatment="Treat",data=X,drop=c("Treat","Post","Post.b"),k2k=F,method="euclidean",cutpoints=list(GTD_SUICIDE=4))
mat
m.data.cem.n <- X[mat$matched,]
sdm.n <- (apply(m.data.cem.n[m.data.cem.n$Treat==1,rightside],2,mean)-apply(m.data.cem.n[m.data.cem.n$Treat==0,rightside],2,mean))/apply(m.data.cem.n[m.data.cem.n$Treat==1,rightside],2,sd)

predata.i <- prematch.data(treatment="SI",lag=12,def=2,type="I",bin=1,data=data)
X <- predata.i[,c("Post","Treat",rightside)]
X <- na.omit(X)
mat <- cem(treatment="Treat",data=X,drop=c("Treat","Post","Post.b"),k2k=F,method="euclidean",cutpoints=list(GTD_SUICIDE=4))
mat
m.data.cem.i <- X[mat$matched,]
sdm.i <- (apply(m.data.cem.i[m.data.cem.i$Treat==1,rightside],2,mean)-apply(m.data.cem.i[m.data.cem.i$Treat==0,rightside],2,mean))/apply(m.data.cem.i[m.data.cem.i$Treat==1,rightside],2,sd)

save(predata.i,predata.n,m.data.cem.i,m.data.cem.n,file="Output/CEM_Data2_RR.RData")


#####################
## Table 7: CEM plot
######################

load("Output/CEM_Data2_RR.RData")

# Model-based estimate
attform <- as.formula("Post ~ Treat + Pre + GTD_SUICIDE + HOLIDAY + POP + CHKINT_NEAR + DIST_MIL + REFUGEE_MIN + DEPORTED")

mod <- glm(attform,data=m.data.cem.i,family="poisson")
pred.y <- predman2(mod=mod,var="Treat",vals=c(0,1))
pred.b <- list(rr=pred.y$rr,rr.lo=pred.y$rr.lo,rr.up=pred.y$rr.up)

mod <- glm(attform,data=m.data.cem.n,family="poisson")
pred.y <- predman2(mod=mod,var="Treat",vals=c(0,1))
pred.c <- list(rr=pred.y$rr,rr.lo=pred.y$rr.lo,rr.up=pred.y$rr.up)

jpeg(file="Output/Table7_CEM.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot3(pred.b=pred.b,pred.c=pred.c,cex=.8)
dev.off()


######################
## Mahalanobis
######################

predata.n <- prematch.data(treatment="SI",lag=12,def=2,type="N",bin=1,data=data)
X <- predata.n[,c("Post","Treat",rightside)]
X <- na.omit(X)
m.out <- matchit(form, data = X,method = "nearest",distance="mahalanobis",discard="both",replace=T,exact=exacts,m.order="random")
m.data.mh.n <- match.data(m.out)
summary(m.out,standardize = TRUE)

predata.i <- prematch.data(treatment="SI",lag=12,def=2,type="I",bin=1,data=data)
#pscore <- glm(form,data=predata.i,family="binomial")
#summary(pscore)
X <- predata.i[,c("Post","Treat",rightside)]
X <- na.omit(X)
m.out <- matchit(form, data = X,method = "nearest",distance="mahalanobis",discard="both",replace=T,exact=exacts,m.order="random")
m.data.mh.i <- match.data(m.out)
summary(m.out,standardize = TRUE)


# Model-based estimate
attform <- as.formula("Post ~ Treat + Pre + GTD_SUICIDE + HOLIDAY + POP + CHKINT_NEAR + DIST_MIL + REFUGEE_MIN + DEPORTED")

mod <- glm(attform,data=m.data.mh.i,family="poisson")
pred.y <- predman2(mod=mod,var="Treat",vals=c(0,1))
pred.b <- list(rr=pred.y$rr,rr.lo=pred.y$rr.lo,rr.up=pred.y$rr.up)

mod <- glm(attform,data=m.data.mh.n,family="poisson")
pred.y <- predman2(mod=mod,var="Treat",vals=c(0,1))
pred.c <- list(rr=pred.y$rr,rr.lo=pred.y$rr.lo,rr.up=pred.y$rr.up)

jpeg(file="Output/Table7_MH.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot3(pred.b=pred.b,pred.c=pred.c,cex=.8)
dev.off()



############################
## Propensity score
############################

predata.n <- prematch.data(treatment="SI",lag=12,def=2,type="N",bin=1,data=data)
X <- predata.n[,c("Post","Treat",rightside)]
X <- na.omit(X)
m.out <- matchit(form, data = X,method = "nearest",distance="logit",discard="both",replace=T,exact=exacts,m.order="random")
m.data.ps.n <- match.data(m.out)
summary(m.out,standardize = TRUE)

predata.i <- prematch.data(treatment="SI",lag=12,def=2,type="I",bin=1,data=data)
X <- predata.i[,c("Post","Treat",rightside)]
X <- na.omit(X)
m.out <- matchit(form, data = X,method = "nearest",distance="logit",discard="both",replace=T,exact=exacts,m.order="random")
m.data.ps.i <- match.data(m.out)
summary(m.out,standardize = TRUE)

# Model-based estimate
attform <- as.formula("Post ~ Treat + Pre + GTD_SUICIDE + HOLIDAY + POP + CHKINT_NEAR + DIST_MIL + REFUGEE_MIN + DEPORTED")

mod <- glm(attform,data=m.data.ps.i,family="poisson")
pred.y <- predman2(mod=mod,var="Treat",vals=c(0,1))
pred.b <- list(rr=pred.y$rr,rr.lo=pred.y$rr.lo,rr.up=pred.y$rr.up)

mod <- glm(attform,data=m.data.ps.n,family="poisson")
pred.y <- predman2(mod=mod,var="Treat",vals=c(0,1))
pred.c <- list(rr=pred.y$rr,rr.lo=pred.y$rr.lo,rr.up=pred.y$rr.up)

jpeg(file="Output/Table7_PS.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot3(pred.b=pred.b,pred.c=pred.c,cex=.8)
dev.off()



############################
## Pre-Matching
############################

predata.n <- prematch.data(treatment="SI",lag=12,def=2,type="N",bin=1,data=data)
predata.i <- prematch.data(treatment="SI",lag=12,def=2,type="I",bin=1,data=data)

# Model-based estimate
attform <- as.formula("Post ~ Treat + Pre + GTD_SUICIDE + HOLIDAY + POP + CHKINT_NEAR + DIST_MIL + REFUGEE_MIN + DEPORTED")

mod <- glm(attform,data=predata.i,family="poisson")
pred.y <- predman2(mod=mod,var="Treat",vals=c(0,1))
pred.b <- list(rr=pred.y$rr,rr.lo=pred.y$rr.lo,rr.up=pred.y$rr.up)

mod <- glm(attform,data=predata.n,family="poisson")
pred.y <- predman2(mod=mod,var="Treat",vals=c(0,1))
pred.c <- list(rr=pred.y$rr,rr.lo=pred.y$rr.lo,rr.up=pred.y$rr.up)

jpeg(file="Output/Table7_PreMatching.jpg",height=.7,width=2.75,units="in",res=150)
rr.plot3(pred.b=pred.b,pred.c=pred.c,cex=.8)
dev.off()



####
# Table 0: Summary stats on matched data
####

match.stat <- as.data.frame(matrix(NA,nrow=4,ncol=8))
row.names(match.stat) <- c("Pre-Matching","Mahalanobis","Propensity Score","CEM")
colnames(match.stat) <- c("Islamist T","Islamist C","Islamist SDM","Islamist % improve","Nationalist T","Nationalist C","Nationalist SDM","Nationalist % improve")
mdatas <- list(list(predata.i,predata.n),list(m.data.mh.i,m.data.mh.n),list(m.data.ps.i,m.data.ps.n),list(m.data.cem.i,m.data.cem.n))
i <- 2
for(i in 1:4){
  sdm.prei <- mean(abs((apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside],2,mean)-apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==0,rightside],2,mean))/apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside],2,sd)))
  match.stat[i,1] <- sum(mdatas[[i]][[1]]$Treat==1)
  match.stat[i,2] <- sum(mdatas[[i]][[1]]$Treat==0)
  match.stat[i,3] <- sdm.prei
  match.stat[i,4] <- 100*(match.stat[1,3]-sdm.prei)/match.stat[1,3]
  
  sdm.pren <- mean(abs((apply(mdatas[[i]][[2]][mdatas[[i]][[2]]$Treat==1,rightside],2,mean)-apply(mdatas[[i]][[2]][mdatas[[i]][[2]]$Treat==0,rightside],2,mean))/apply(mdatas[[i]][[2]][mdatas[[i]][[2]]$Treat==1,rightside],2,sd)))
  match.stat[i,5] <- sum(mdatas[[i]][[2]]$Treat==1)
  match.stat[i,6] <- sum(mdatas[[i]][[2]]$Treat==0)
  match.stat[i,7] <- sdm.pren
  match.stat[i,8] <- 100*(match.stat[1,7]-sdm.pren)/match.stat[1,7]

}
match.stat


write.table(round(match.stat,2),file="Output/Table0_MatchStat.txt")
write.table(print(xtable(match.stat)),file="Output/Table0_MatchStat.tex")




####
# Table 0: Detailed summary stats on matched data
####

rightside <- c("Pre","PreG","POP","REFUGEE_MIN","DIST_MIL","DEPORTED","ELEVATION","SLOPE","LANGUAGE","FOREST","CHKINT_NEAR","GTD_SUICIDE","HOLIDAY","LONG","LAT","YEAR","MONTH","WID")
length(rightside)
labz <- c("Insurgent violence (pre-T)","Government violence (pre-T)","Population density","Distance to nearest refugee camp","Distance to military base","Deported in 1944","Elevation","Slope","Percent Russian speaking","Forest","Distance to international checkpoint","Global suicide terrorism","Muslim holiday","Longitude","Latitude","Year","Month","Week")
cbind(rightside,labz)
apply(predata.i[,rightside],2,function(x){mean(x,na.rm=T)})

mdatas <- list(list(predata.i,predata.n),list(m.data.mh.i,m.data.mh.n),list(m.data.ps.i,m.data.ps.n),list(m.data.cem.i,m.data.cem.n))

i <- 1
mchtz <- c("PreMatching","Mahalanobis","PropensityScore","CEM")
for(i in 1:4){

  match.stat.i <- as.data.frame(matrix(NA,nrow=length(rightside),ncol=5))
  row.names(match.stat.i) <- labz
  colnames(match.stat.i) <- c("Treated","Control","SDM","TTest","KSTest")
  match.stat.n <- as.data.frame(matrix(NA,nrow=length(rightside),ncol=5))
  row.names(match.stat.n) <- labz
  colnames(match.stat.n) <- c("Treated","Control","SDM","TTest","KSTest")
  
  match.stat.i[,1] <- apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside],2,mean)
  match.stat.i[,2] <- apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==0,rightside],2,mean)
  match.stat.i[,3] <- abs((apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside],2,mean)-apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==0,rightside],2,mean))/apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside],2,sd))
  for(j in 1:length(rightside)){match.stat.i[j,4] <- t.test(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside[j]],mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==0,rightside[j]])$p.value}
  for(j in 1:length(rightside)){match.stat.i[j,5] <- ks.test(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside[j]],mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==0,rightside[j]])$p.value}
  
  match.stat.n[,1] <- apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside],2,mean)
  match.stat.n[,2] <- apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==0,rightside],2,mean)
  match.stat.n[,3] <- abs((apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside],2,mean)-apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==0,rightside],2,mean))/apply(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside],2,sd))
  for(j in 1:length(rightside)){match.stat.n[j,4] <- t.test(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside[j]],mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==0,rightside[j]])$p.value}
  for(j in 1:length(rightside)){match.stat.n[j,5] <- ks.test(mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==1,rightside[j]],mdatas[[i]][[1]][mdatas[[i]][[1]]$Treat==0,rightside[j]])$p.value}

  write.table(round(match.stat.i,2),file=paste("Output/Table0_DetMatchStat_",mchtz[i],"_I.txt",sep=""))
  write.table(print(xtable(match.stat.i,digits=3)),file=paste("Output/Table0_DetMatchStat_",mchtz[i],"_I.tex",sep=""))

  write.table(round(match.stat.n,2),file=paste("Output/Table0_DetMatchStat_",mchtz[i],"_N.txt",sep=""))
  write.table(print(xtable(match.stat.n,digits=3)),file=paste("Output/Table0_DetMatchStat_",mchtz[i],"_N.tex",sep=""))

}



########################
## Poisson tables
########################

attform <- as.formula("Post ~ Treat + Pre + GTD_SUICIDE + HOLIDAY + POP + CHKINT_NEAR + DIST_MIL + REFUGEE_MIN + DEPORTED")
mdatas <- list(list(predata.i,predata.n),list(m.data.mh.i,m.data.mh.n),list(m.data.ps.i,m.data.ps.n),list(m.data.cem.i,m.data.cem.n))

mod1 <- glm(attform,data=predata.i,family="poisson")
mod2 <- glm(attform,data=predata.n,family="poisson")
mod3 <- glm(attform,data=m.data.mh.i,family="poisson",weights=m.data.mh.i$weights)
mod4 <- glm(attform,data=m.data.mh.n,family="poisson",weights=m.data.mh.n$weights)
mod5 <- glm(attform,data=m.data.ps.i,family="poisson",weights=m.data.ps.i$weights)
mod6 <- glm(attform,data=m.data.ps.n,family="poisson",weights=m.data.ps.n$weights)
mod7 <- glm(attform,data=m.data.cem.i,family="poisson")
mod8 <- glm(attform,data=m.data.cem.n,family="poisson")



varlabels <- c("Selective tactics (T)","Violence (pre-T)","Global suicide terrorism","Muslim holiday","Population density","Distance to border crossing","Distance to military base","Distance to nearest refugee camp","Deported in 1944")
stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,style="apsr",covariate.labels=varlabels)





